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Abstract 


Theoretical models are developed and numerical studies are conducted on various 
types of flows including both elliptic and parabolic flows. The purpose of this study 
is to find better higher order closure models for the computations of complex flows. 

This report summarizes three new achievements: i) Completion of the Reynolds- 
stress closure by developing a new pressure-strain correlation, ii) development of a 
parabolic code to compute jets and wakes, and iii) application to a flow through 180 ° 
turn around duct by adopting a boundary fitted coordinate system. 

In the above mentioned models near-wall models are developed for the pressure- 
strain correlation model and third-moment, and incorporated into the transport equa- 
tions. This addition improved the residts considerably and is recommended for future 
computations. 

A new parabolic code to solve shear flows without coordinate transformations is 
developed and incorporated in this study. This code uses the structure of the finite 
volume method to solve the governing equations implicitly. The code was validated 
with the experimental results available in literature. 
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constants used in the source terms of the turbulent energy 
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constant used in of < UiUjUk > 

wall constant used in <f>ijk,T of < UiUjUk > 

constant used in calculation of turbulent eddy viscosity 

coefficient for the dissipation rate of < UiUjUk > 
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Constants used in algebraic equations of < Ui 0 > of Launder and 
Samar aweera 

constant used in P ijkx of < UiUjU k > 
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constants used in boundary conditions of < UiUj > 

skin friction coefficient, = t- t v 

2^ in 

constant used in the near wall model of < UiUju k >, = kC^ 
wall static pressure coefficient, = Zr Pr *t 
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constant pressure heat capacity 
constant used in diffusion rate of < U{Uj > 

Temporary constants used in term IV of < mujO > 
nozzle diameter 
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diffusion rate of < urf > 
diffusion rate of < muj > 
diffusion rate of < um : 0 > 
diffusion rate of < UiUjU k > 

height of channel downstream of the step, = [Y 0 + H) 

generation rate of turbulent kinetic energy 

primary generation or production rate of < Ui Uj > 

step height of backward-facing step 

secondary generation or production rate of < U{Uj > 

turbulence intensity 

turbulent kinetic energy 

mixing length 

pressure fluctuation 

mean pressure 

cell Peclet number 

production rate of < U{6 > 

generation or production rate of < u iUj 0 > due to mean strain rate 
generation or production rate of < u iUj 0 > due to interactions 
between Reynolds stresses and second moments of velocity- 
temperature products 

generation or production rate of < > due to mean strain rate 

generation or production rate of < Ui u 3 u k > due to interactions 
between the Reynolds stresses and their gradients 
Prandtl number 

turbulent Prandtl number 
wall heat flux per unit width 
radius of the nozzle 
orifice Reynolds number 
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step-height Reynolds number, =^R 

source terms of transport equations in the discretization equation 
— Su -f- Spcf>p 

linearized source terms in S $ 
mean temperature 
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inlet stream velocity 

Maximum temperature in shear layer 

local wall temperature 

fluctuating velocity in x-direction 

second moments of temperature- velocity products 

axial component of < U{Uj > 

Reynolds shear stress 
Reynolds Stresses 

third moments of velocity-temperature products 

hydrodynamics third moments of turbulence 

mean velocity in x-direction 

friction velocity 

mean velocity components 

discharge velocity in coflowing jets 

inlet velocity 

Maximum velocity in shear layer 
surrounding velocity in coflowing jets 
fluctuating velocity in y-direction 
transverse component of < U{Uj > 
mean velocity in y-direction 
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normalized coordinate, — 
normal distance from the wall 
reattachment length measured from the step 
half-maximum temperature width 
half-maximum velocity width 
coordinate normal to the wall 
coordinate parallel to the wall 
channel width upstream of step 
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kronecker delta 

dissipation rate of turbulent kinetic energy 
dissipation rate of < > 

dissipation rate of < U{Uj > 
dissipation rate of < UiUjO > 
dissipation rate of < UiUjUk > 

similarity coordinate for temperature in shear layer, = ^ 
similarity coordinate for velocity in shear layer, = ^ 
diffusion coefficient in general governing and discretization 
equations of parabolic equation 

diffusion coefficients in x- and y-direction in general governing and 
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effective kinematic viscosity, = v + v t 
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dependent variables in general governing and discretization equat.ons 
pressure-heat flux effects of < u$ > 

wall correction of pressure-heat flux effects of < utf > 
pressure-strain correlations of < UiUj > 

pressure-strain correlations of < > due to fluctuating components 

pressure-heat flux effects of < U{UjO > 
pressure-strain correlations of < UiUjUk > 
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SUBSCRIPTS 
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IN 

inlet station conditions 
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node cell P 
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reattachment station conditions 
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cell boundaries at east, west, north and south of cell P 
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reference station (x/H = -4) conditions 
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Chapter 1 


Introduction 


In the past decade, Computational Fluid Dynamics (CFD) is assuming an im- 
portant role in design with the advancements in computer technology and numerical 
method. Peterson and Bailey [1] cited that supercomputers are becoming so powerful 
that aerodynamic flow simulations are almost as good as wind tunnel testing. It was 
noted that the sophistication of computational fluid dynamics is so advanced that 
results of computed surface pressure on the shuttle launch configuration with Mach 
number 1.55 and Reynolds number of 2.5 x 10 6 coincided with a wind tunnel testing 
and flight data. 

Most commercial firms and research institutes use the conventional k — e model of 
turbulence which gives satisfactory results for planar, wall-shear layer and favorable 
pressure gradient internal flows. This model uses the eddy viscosity concept assuming 
that turbulence is locally isotropic or the normal components of the Reynolds stresses 
have equal magnitudes in all directions. However, in flow fields where large rates of 
strain are imposed by curvature, sharp corner, swirl and body forces, the magnitudes 
of Reynolds stresses are no longer equal in all directions. Thus, the isotropic assump- 
tion of the k — e model fails to give realistic results. Driver and Seegmiller [2] pointed 
out that the k — e model of Jones and Launder [3] overestimated the isotropic tur- 
bulent viscosity in recirculating flows resulting a higher spreading rate for the shear 
layer giving a premature reattachment. 

The next attractive alternative to correct this problem is to use the algebraic 
stress models (ASM) which are approximations of the full transport equations of the 



Reynolds stresses in algebraic form. Hutchings and Iannuzzelli [4] reported compu- 
tations of highly swirling flow with the k — £ model and the algebraic stress model 
and illustrated the deficiency of the k — e model in the predictions of swirling flows. 
The disadvantage of using the algebraic stress model is that it does not take the 
convective and diffusive transport of Reynolds stresses into considerations. In high 
Reynolds number flows, the convective processes play an important role and should 
not be omitted in the computations. 


1.1 Literature Survey 

Chandrsuda and Bradshaw [5] carried out the experiment with a step-height 
Reynolds number of about 10 5 , a step height of 51 mm, inlet nozzle of 127 mm and 
inlet freestream velocity of 31.5 m/s . The reattachment point was found from surface 
oil-flow measurements to be about 5.9 step heights downstream of the step. Hot-wire 
and pressure-probe were used to measure the important quantities, while the skin- 
friction coefficient was obtained from surface tube measurements. Although there 
were no specific number given for the uncertainties of measurements, measurements 
inside the recirculation zone were recommended for quantitative use only. 

Driver and Seegmiller [6] conducted the experiment with a step-height Reynolds 
number of 3.78 x 10 4 , a step-height of 1.27 cm, inlet nozzle of 10.16 cm and freestream 
inlet velocity of 44.2 m/s. The reattachment point was found to be 6.2 step heights 
downstream of the step. The uncertainty of the measured wall static pressure coeffi- 
cient and skin-friction coefficient were =b 0.0009 and ± 8% (± 15 % in the separated 
region of the flow) respectively. No uncertainty values were reported for other quan- 
tities. 

Eaton and Vogel [7] performed the experiment with a constant heat-flux surface 
behind a single-sided backward-facing step . It was conducted with a step-height 
Reynolds number of 28000, a step height of 3.8 cm, freestream inlet velocity of 11.3 
m/s and constant heat flux of 270 W/m. The uncertainties of the heat flux, mean ve- 
locity, and fluctuating velocity were 1% , 1% and 2%, respectively. The reattachment 
point was 6| step heights downstream of the step. 

Figure 1 shows the experimental conditions for all three experiments. 
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(a) Chandrsuda and Bradshaw 

H =57 mm 
Yo = 127 mm 
Uin = 31.5 m/s 
Re H = 10 s 
x R = 5.9 H 

(b) Driver and Seegmiller 

H = 1.27 cm 
Yo — 10.16 cm 
U in = 44.2 m/s 
Re J{ = 3.78 x 10 4 
XR = 6.2 H 

(c) Eaton and Vogel 

H = 3.8 cm 
Yo =15 cm 
R in = 11.3 m/s 
(q)l = 270 W/m 
Rch = 28000 
x R = 6| H 


Figure 1: Experimental conditions for all three backward-facing step 
investigated 


experiments 



Antonia et al. [8, 9, 10, 11, 12, 13] conducted the experiment with a variable- 
speed centrifugal squirrel cage blower which supplies air to a setting chamber followed 
by a vertical nozzle of contraction ratio 20:1. The nozzle exit velocity was 9 m/ 5 , 
orifice Reynolds number, Red, was 7550 and jet temperature was maintained at a 
nominal temperature of 25°C' relative to the ambient temperature. Measurements 
for mean velocity and temperature were done using a 5 fim hot wire with a DISA 
55M01 constant temperature anemometer and a 0.63 fxm cold wire drawing 0.1 mA in 
a constant current circuit. Fluctuating temperatures were measured with cold wires 
operated at very low overheat as resistance thermometers in constant current circuits. 

Heskestad [14] performed the experiment by exhausting a constant velocity jet 
into still air. The exit nozzle velocity and orifice Reynolds number were 1.27 cm 
and 3.4 x 10 4 , respectively. Hot wire was used to measure all quantities and true 
self-preservation was not attained [15]. 

Gutmark and Wygnanski [16] measured mean velocities, turbulence intensities, 
third- and fourth-moment, as well as, two-point correlations and the intermittency 
factor by hot wire. Measurements were made up to a distance of x/d = 120 and 
flow was found to be self-preserving beyond x/d = 40. The exhaust velocity, orifice 
Reynolds number and exit nozzle width were 35m/ 5 , 3 x 10 4 and 1.3 cm, respectively. 

Everitt and Robins [15] performed an experiment on submerged jet with a variable 
nozzle width ranging from 0.32 cm to 2.54 cm. Documentations were made on the 
Reynolds stresses and triple-moment products of velocities. 

Wygnanski and Fiedler [17, 18] found that most previous investigators did not 
measure the turbulent quantities far downstream enough to ensure self-preservation 
was achieved. Measurements were made beyond 100 nozzle diameters by the authors 
to ensure self-preservation conditions were attained 

Numerical studies of separating and reattaching flows are less numerous than 
experimental studies. During the sixties, Gosman et al. [19] pioneered the theoretical 
work on turbulent flows predictions. Pope and Whitelaw [20] studied the near wake 
flows using three different models of turbulence. 

Schlichting [21] attacked the problem of plane laminar jet and computed the 
velocity across the jet by an approximate numerical method from the fundamental 
equations of constant density viscous flow. In going through Schlichting’s work, Bick- 
ley [22] and Goldstein [23] found that in a plane jet case, the equations are integrable 



in close form. 

In the past few years, increase attentions have been given to the development 
of second-order turbulence models. This is the simplest closure level which can 
incorporate the essential turbulent flow characteristics such as transport, pressure- 
interactions, dissipations and effects of external force fields. 

Hanjalic and Launder [24] performed computational works on a plane mixing 
layer, plane jet, boundary layer and channel flow by using second-moment turbulence 
closure. Dekeyser [25] computed an asymmetrically heated coflowing plane jet with 
second-moment closure. Sini and Dekeyser [26, 27] performed numerical work on 
turbulent shear flows using second-moment turbulence closure and on turbulent plane 
jets and forced plumes using the k - e model of turbulence, respectively. Dekeyser 
and Launder [28] modeled the triple moments of velocity and velocity-temperature 
of a coflowing jet. Samaraweera [29] documented numerical results of two and three 
dimensional temperature field and later combined with Launder [30] to apply second- 
moment turbulence closure to investigate heat and mass transport in thin shear flows. 

All computations were performed by using the computer code of Patankar and 
Spalding [31] with some minor modifications of staggering the < uv > cell. 


1.2 Objectives 

The objectives of this study are listed below : 

1. Refine the existing Reynolds stresses and third-moment turbulent computations. 

2. Develop a new parabolic code for the computations of parabolic flows without 
coordinate transformations 

3. Apply the above formulations to a 180° turn around duct flows by using a 
boundary fitted coordinates. 
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Chapter 2 

Mathematical Formulations 

This chapter defines all the transport and algebraic equations governing the de- 
pendent variables used in this study. 


2.1 Transport Equations 

Reynolds-Stress Equations 


u d<uu J > = _ + ^ + Di , 

OXk 


where 


dUi 


G ‘i = -{ <u ‘ ut> dT k +<UiUk> dx k 


duA 

dx k ) 


D ” - c 4 

Cxj — 


k d < UiUj > 

- < U k Ul > r 

e axi 


There are a few models for <j>. 




Model 1 : Naot et al. [34]: 


( 1 ) 


<j>ij = -C + 2 ( Gij - + <j>ij , i 


(2) 
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Model 2 : Amano et al:(see section 2.3) 


- (7C* - 10)G Dt(g + @ 

-2(C* - 1) (a„ - |m?) - 2 (C* - 1) (h v - |«»g) + (3) 


Model 3 : Launder et al: [36]: 


4*ij — 


(<?* + 8) 


(*-! 


SaG ) - 


(8Q 2 - 2) 


fiii - 


where 


(30^-2) /at/, dU 3 

55 \dxj dxi 


TT { dU k dU k 

Hij = - I < UiU k > ^-+ < u i u k > 0^7 

~ / dUi dUj 

hij — I UjUk “1“ UiUk 


$ ij ,\ ^<^1 ^ ^ ^ ^ 

„ at/< 

G = - < UiUj > - — 


Third-Moment Equations : Low- Reynolds number model [37]: 

a (< Ill'll jtlk > ) — Pijki " 4 * Ptjkj H" (pijk "I" $ijk,w &tjk Dijk 

OX l v v / 


where 

D ~ / at/* , at/,- at/, A 

Pijki = -G g I < UjUjU; > < UjUkUl > + < «Jt «,«! > -7^ I 

0 / a < u,u, > a < UjUk > a < > 

Pijk 2 = - I ■ < WfcW/ > 1 — + < UiU, > + < UjUl > _^_2 — 

M ^ < «,•«>«* > „ Art n /aA:^ 2 ,, 

^tjkyT - c \ fc l|€ ’ c ^cw 2l/ ay J 11 

& ijk — C^sk l 


Diik = 


a / a 


J/— < U.-UjU* > 


dxi\dxi ,]l 
High- Reynolds number model: 


< U.UjUfc > 

<Pijk,T = ~ C 7 < 

£ijfc = 0.0 
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Heat Transfer Governing Equations 


2.2.1 Scalar Double Product Equations 


The transport equations of second moments of velocity-temperature, < u,0 >, are 
formulated and simplified through closures. After neglecting insignificant terms, the 
transport equations reduce to the following [38]: 


d_ 

'dxj 


Ui-£r{<»iO >) 


dT . dUi 

< u ‘ u > > < u ’ 6 > 


d 


dxi 


. .d <UiO > a 

(a + v) 5- < UiUjO > 


dxj 


ii 


pdO , , x ^ dd dui 

+ < -7T- > -(« + v) < ^ — n — > 

p dxi dxj dxj 


(7) 


in 


IV 


which can be written as 

U j — (< Uid >) = Pie + Die + 4>ie + — £»« (®) 

dxj 

where 

P ie = Production rate of < mO >. (Term I) 

Dig = Diffusion rate of < U{0 >. (Term II) 

(f> i9 = Pressure-heat flux effects. (Term III) 

<f>ie,-w — Wall correction of pressure-heat flux. (Term III) 

Eie = Dissipation of < >. (Term IV) 

Term I is explicit in character and thus needs no further modifications. The 
diffusive rates, which is term II contains the triple products < UiUj0 > which has to 
be closed by using an appropriate closure, is decomposed into simpler form following 
Launder’s procedure [39]. 


k. d<U{ 0 > 

< UiUjO >= -0.15(2-) < UjU, > ^ 

The pressure- heat flux effects can be modeled as follow [39]. 

dUi 


(9) 
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The near-wall correction of pressure-heat flux effect is proposed by Launder and 
Samaraweera [30] as 


8 } w — 


-0.1- < UiO > -0.02 < ui9 > 

k 


(# - B 

y (JXi C/Xi J 


k §_ 

£X n 


( 11 ) 


The dissipation rate is assumed to be negligible in accordance with the assump- 
tions made in the Reynolds stresses closures. 

The final form of the transport equations of < U{6 > is 


U j — ( < u t 0 > ) = Pie + Die + <t>ie + <j>ie,w 

axj 


( 12 ) 


Pie 

Die = 
<f>ie = 
4^i0 t w 


( dT , dUi\ 


d_ 

dxj 


d < > 0.3A; d < U{0 > 

(a + v) — b — < UjUi > 


dxj 


dxi 


dU' 

-CiO 1 — < UiO > -\-CiO ,2 < u l@ > 

’ k oxi 


—0.1— < u 

k 


„ (dUi dUA] k § 
,0 > -0.02 < u,e > (4^— - J — 


2.2.2 Scalar Triple Product Equations 

In order to evaluate the triple products of velocity-temperature fluctuations, 
< UiUjO >, transport equations were formulated and are given as follows : 


£/^— (<C UiUjO >) 

ox k 


< UiUjUk > 


dT 

dx k 


n dUi Q dUA 

< u ‘ Uk6 > dI~T < UkU ‘ > ) 


() & 

[< UiUj— — {u k 0— < u k 0 >) > + < Ui0—{ujU k - < UjUk >) > 

dx k VX k _ 


II 


d 


+ < Uj0- - — ( U k Ui - < U k Ui >) >] 

dx k 


ii 


f d , dO x n d , fi duj s ^ n d r fidui^J 

< > + < Ui W7fe' + W > , 


III 
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' U{9 dp u 

< — — £~ > + < 

K P OXj p 


£ dp \ 

p dxi j 


(13) 


IV 


which can be written as 

0 

U k -z — (< UiUjO >) = Pije , i + Pijea + Dije + - eije (14) 

OXk 

I = Production rate due to the mean strain rate and temperature 

gradients 

II = Production rate due to the interactions of < U{Uj > and 

gradients of < UiUj > with the heat flux components, < U{9 > 

III = Diffusive and dissipative effects due to molecular viscosity 

IV = Pressure-heat flux effects 

In closing the above equations, term I needs no further modifications since it is 
explicit. Term II can be rearranged and written as: 

d < u k 6 > a^ d< u i Uk > , / a ^ 9 < UkUi > 

II — < U{Uj > h < Uid > — h < UjV > 


_d_ 

dxi 


dx k 

(< UiUjU k 0 >) 


dx k 


dxi 


(15) 


The quadruple terms are assumed to be Gaussian and can be split up as 

< UiUjU k 0 > — < upij >< u k 0 > + < UjUk >< U{0 > + < u k Ui >< Uj9 > (16) 

Differentiating equation ( 16 ) with respect to x k and substituting the result into 
equation ( 15 ) yields 


TT d < Ui9 > d < Ujd > . d < UiUj > 

II — — < UjU k > £ < u k Ui > 5 < u k 6 > 


dx k K ’ dx k 

Term III can be rearranged and cast into the following form : 


dx k 


(17) 


III = (a + 2v) 
-a (2 < 


d ( d 


dx k \dx k 
d0 duiUj 
dx k dx k 


< UiUjO 


> + <0 


( dui dujO 
-v[2< — ~—£— > + <« 
y dx k dx k 

f duj duiO d 2 UiO* 

-H 2 < — ~ ” > + < Uj ^ „ > 


dxi > ) 

£^£1 A 

• dxi > ) 


dx k dx k 


dx\ 


( 18 ) 
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The first term in equation ( 18 ) represents the laminar diffusion while the rest of 
the terms express the dissipative effects. Terms with asterisks (*) consist of second 
derivatives of the second moments and are assumed to be negligible. Empirical coef- 
ficients were optimized to compensate the effects of the asterisk terms. Term IV is 
the pressure-heat flux effects term, and can be modeled as [40]: 


IV « 


^ UjU^O ^ 


+- < U k Uit 



CVy < U,UjO 



(19) 


where the second term in equation ( 19 ) is merged into term I of equation ( 13 ) and 
the coefficient C ei was adjusted accordingly. 

The final form of the transport equations are given below as: 

d 


where 

Pij8,\ 
PijB, 2 
D tJ e 

£ij& 


Uk 


dxi 


(< UiUjO >) — PijB, 1 + PijB , 2 + Dijo + — £ijg 


(20) 


dT „ dUi A dUj 

< UiUjUk > 7J h < UjUkO > — 1- < UkUiV > 


< UjUk > 


dx k 

d <Ui& > 


dxi 


dx k 


dxi 


f < U k Ui > 


d < Uj8 > . d < UiUj > 

3 + <u k 8> 3 


dxi 


dxi 


d 

dx k 


d 


(o- + 2is )~ — < UiUjO > 
dxk 


—Cg-iT < UiUjO > 
n e n 

Cg £ j < U k 6 > -r — 
k oxk 


2.3 Pressure Strain Correlations 


2.3.1 Formulation of a Pressure Strain Correlation 

Following Chou [41], a Poisson equation for the fluctuating pressure, p, can be ob- 
tained by taking the divergence of the equation for the turbulent fluctuating velocity, 
u,-, which can be re-expressed to write the pressure-strain correlation as: 


p dui 
P 


>_ Sii + s 



dvol 

~R~ 


( 21 ) 
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where 



( 22 ) 

(23) 

(24) 


where the first and second terms in the integral are denoted as <f>ij i i and <f>ij t 2 , respec- 
tively, and Sij is a surface integral which is negligible away from a solid wall. Sij term 
is also developed and described in the next subsection. 

Following Rotta’s proposal [42], i was formulated same as one given by equation 
( 5 ) with Cfa = 1.5. Furthermore, may be approximated as 

fo = (25) 

where 11™* is a fourth-order tensor which satisfies the kinematic constraints of 


TT77U TTtm TT 1771 

- AA /j “ ll jl 


(26) 


nr = 0 


(27) 


n^‘ = 2 < u m ui > 

(28) 

njj = 2 < it[Uj > 

(29) 


Equation ( 25 ) suggests that II jj* can be approximated by a linear combination of 
the Reynolds stresses as: 

njj * = Aj < U m Ui >< UlUj > +B^r < U m Ui >< UiUj > 

nl AC 

+c \ < HmHj ^ ^ U%Ul ^ <C Ufjilti 

+ ^ 11 {'ll j ^ U {11 1 

^ Hm,Hj ^ ^ H mil l y> ) 

^ Ill'll j ^ H" “I” 


(30) 
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where the condition of equation ( 26 ) is already imposed. The application of equations 
( 27 ) through ( 29 ) enables two of these constants to be expressed in terms of the 
third, namely Cj> 2 and can be written as 


A = B = —C = 


10 - IC^ 
2 


p \dxj dxi 


D = C * 2 
E = C+,-2 

L __ _ (i±fZfe) 
m - ( 5 = 32 .) 

The final form is given as 

p ( du{ duj\ t ( 2 \ 

< ~p{d^ + d^) > ~ ~ C *k UiUi > ~3 6ij V 

+ (7C* - 10)G - §*«) 

- 2 (C* - 1)(G* - |%G) 

- 2(C*-1 )(//,-,- 

- < 31 » 

The coefficient 0 $ 2 can be determined by simplifying the Reynolds stresses described 
by equation ( 1 ) into an algebraic equation which was originally developed by Rodi 
[43]. 

If the convection-diffusion string of the < upij > is set equal to that of the 
turbulence kinetic energy, k , it follows that, 


q — (Ui < UiUj >) - Dij 


< U{Uj > [ d 


-( u,k ) - D k 


where D,j and I)' K denote, respectively, the diffusion rates of < U{Uj > and k. After 
some manipulation, equation ( 32 ) can be re-arranged into the following form, 


< u^j > 2 

— i s Sii = 


(2 - 3 )(G, V - ISjjG) + (2 Cy, - 2)(g„ - j frg) 
e -f- 10G — 7(7^ G 
ACfo — 2 \ & /dt/, 

4- 10G — 7C<f, 2 G J 5 Vdx, + dxi) 
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By assuming that the dissipative action is isotropic, and the flow is parabolic, equation 
( 33 ) can be further reduced to : 


< uu > 2 _ (4(7*2 — 8) 

k 3 3(C^ 1 — 7 (7* 2 4- 10) 


= 0.3 


(34) 


< vv > 2 _ (4C* 2 — 2) 

k 3 3((7*, — 7C* 2 + 10) 


= -0.18 


(35) 


< ww > 2 _ (10 — 8(7* 2 ) 

k 3 = _ 3(C'^ 1 - 7C* 2 + 10) 


= - 0.12 


(36) 


The above set of relations are compared with typical experimental data for several 
parabolic flows, and the coefficient C^ 2 is found as: 

C* 2 = 1.05 

2.3.2 Formulation of a Near Wall Model 

In parallel to the former derivation for the pressure-strain correlation the near-wall 
correlation can be set as: 




c;|<< u,u, > -§««*) + + nri 


/(- 

y 


(37) 


where the function f{l/y) represents the wall controlling function which enforces 
the total term become negligible in the region away from a solid wall. If the forth- 
order tensor f is expanded, we have 


nr 


= A'^r < U m Ui >< UlUj > +5 — < U m Ul >< UiUj > 

(7 — ^ U{Ul D U m Ui 

K 

+ E (& m l < U^IL j > <C > 

T &H tl m Uj > ^ V’mV'l ^ ) 

+ Cfa8 m i < U{Uj > + L SraiSljk + 


(38) 



where the following conditions are imposed. 


nr = 0 


(39) 


< uu > 2 

k 3 


0.51 


(40) 


< vv > 2 

k 3 


-0.42 


(41) 


< ww > 2 

k 3 


-0.09 


(42) 


< uv > 
k 


= 0.24 


(43) 


In equations ( 39 ) through ( 43 ), these experimental data are introduced from 
Champagne et al. [44] in order to solve the coefficients in terms of C^ 2 . 

Experimental data of Champagne et al. [44] were used to obtain the coefficients 
in equations ( 39 ) through ( 43 ) in terms of Cfo. 


A' = -C' 

B' = 0 
E' = 0.11 Cfr 
D' = — 1.70C* 2 
L' = M' = —0.04C^ 2 


Then the final form of the second term in the parentheses of equation (21 ) becomes: 


dUi 

f)x 


(nr + nr ) - 1 -55C 02 Gij - 0.216^2 G - QAlC^Hij 


—C'foHij - O.OSC^ 


(dlE dUA 

\ dxj dxi J 


(44) 



As before if equation ( 27 ) is substituted into the algebraic stresss equation 
(equation ( 32 )) then, after considerable manipulation. We obtain as: 


<uu> 2 _ (2.67+ 1.56 C^) 

k 3 _ (Cfa — 7C^ 2 + 10 — C'tJ 


= 0.51 


(45) 


< vv > 2 _ (0.67 - 1.76(7* - 2 C*) 

k 3"(Q 1 -7^ 2 + 10-C; i ) 


= -0.42 


(46) 


<ww> 2 _ (2.45(7* - 3.33) 

k 3 " (^-7^ + 10-^) 


= -0.09 


(47) 


< uv > 
k 


0.88(7* - 0.4 


1 0.5 


lC* - 7 <7* + io - c' d 


'<$>2 


4 > iJ 


(3 - 0.45C fe )^ - (C;, - 2 + 2.1 1 


t 0.5 


c t , - vu., + io - c;, 


(48) 


As before by comparing with the experimental data of Champagne et al., coeffi- 
cients are obtained as follow: 


Ci = -4.28 (49) 

C* = 1.18 (50) 

2.4 Constants Used in the Computations 

Hydrodynamics Equations 

Table 1 lists the values of all constants used in the hydrodynamics equations. 
Heat Transfer Equations 


Table 2 lists the numerical values of constants used in heat transfer equations. 


Variable 

Numerical values 

c. 

0.30 

C* 

1.5 

Cfo 

0.4 * 

c v 

0.09 

a, 

3.0 f 

C^yj 

8.0 

C e . y 

0.30 

Cl 

1.44 

C2 

1.92 

Cn 

1.21 

C\2 

-0.24 

C 22 

0.24 

Cl 

2.55 

C a 

0.25 

Oil 

-8.14 x 10" 3 

&2 

-1.72 x 10~ 2 

<*3 

-4.80 x 10- 2 

C*4 

-1.02 x 10 _1 

K 

0.42 


1.0 


* The values of C ^ for parabolic code are: 

= 0.75 for model of Amano et al 
= 0.6 for model of Naot et al 

f The value of C 1 for high- Reynolds number model = 5.8 


Table 1 : Values of Hydrodynamics constants. 


Variable 

Numerical values 

C\T 

0.31 

C2T 

0.16 

Ce-y 

6.0 

Coc 

0.10 

c m 

3.2 

C ie , 2 

0.5 


Table 2: Values of heat transfer constants. 





2.5 Boundary Conditions 


The numerical solutions of the transport equations require the provision of good 
boundary values and this section describes the treatments of boundary nodes used in 
this study. 


2.5.1 Backward- Facing Step Flow 

Inlet Conditions 

The U- velocity is calculated from the step-height Reynolds number and assumed 
to be uniform throughout the inlet section of the computation domain. V-velocity is 
set to zero. 

The turbulent kinetic energy and turbulent energy dissipations are computed from 
the inlet turbulent intensities by the following relations: 


Jr* 

II 

«s> . 

(51) 

k §• 

(52) 

£ = A H 


The inlet Reynolds stresses are given the following values: 



2 

(53) 

< uu > — 

3* 

< vv > = 

3 

(54) 

< uv > = 

— O.OlJfc 

(55) 


The inlet values of the third moments velocity products are specified by the alge- 
braic model of Daly and Harlow [45] as 


_ nr k d < UiUj > 

< UiUjUk >= —0.25— < UkUi > — - — 

J tr s ! nr * , 


(56) 


The inlet temperature is set to the value used in the experiment. While the second- 
moment of velocity-temperature products are specified by assuming that fluctuating 
quantities are 5% of the mean quantities. This can be written as: 


< ufi >= (0.05) 2 T/^f//Ar 


(57) 
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The triple products of velocity-temperature correlations are specified by using 
model of Launder [39] as: 

/> n i ( <9 < UjO > 

< UiUjO >— — 0 . 10 — < U{Ui > 1 - < iijUi > 

£ \ OXi 

Initial Values 

The initial values for U-velocities are set by using the conservation of mass, which 
reduces to UjnYq = U D c for incompressible flows. The V velocities are given zero 
values. 

The turbulent kinetic energy, turbulent energy dissipation and Reynolds stresses 
are given constant values of equal magnitudes with the respective inlet conditions. 

The initial values of triple products of velocity are prescribed by the algebraic 
equation of Daly and Harlow [45]. 

All initial values of heat transfer variables are the same as their respective inlet 
values. 

Wall Boundary Conditions 

At the wall boundary, the wall function treatments are employed for the momentum 
and the turbulent kinetic energy. For the mean momentum equations, the velocity 
gradient is used to determine the wall shear stress which is then used to prescribe the 
boundary values. For the turbulent kinetic energy, the wall shear stress is incorporated 
into the generation rate and thus introducing the wall effects. 

The dissipation rate is evaluated under the local equilibrium condition in terms 
of the turbulent kinetic energy as 


d < UiO > 
dxi 


(58) 


€ = 



Cix n 


(59) 


Launder et al. [36] obtained a correlation between < u l Uj > and U T for channel 
flows as 


< UiUj >= CijU * - (1. - 




Vi dP_ 
P dy 2 


(60) 



20 


The coefficients C,j are given as 

Cm = 5.1 
C \2 = - 1.0 
c 22 = 1.0 

Using the correlation of k and U T in the wall proximity derived by Hanjalic and 
Launder [24]: 

k « 3.5U t 2 (61) 

the boundary values of < iqu,- > can be expressed in terms of k as: 

< U{Uj >= Cijk — (1. — Sij) — — — (62) 

P dy 2 

where the coefficients Cij are now 

Cm = 1-214 
Cm = -0.24 
C 22 = 0.24 

The algebraic equation of Shir [46] was combined with the wall values of < > 

to prescribe the wall values of < U{UjUk > and is given as follow: 

y A P 

< UiUjUk >= 0.04 — — Cijk - (1. - Sij) — — (63) 

J e dx k [ p dy 2 \ 

For the temperature equation, a constant heat flux of 270 w/m is applied along the 

wall downstream of the step. The heat flux is introduced into the solution domain by 

supplementing the source term at the wall adjacent cell for the temperature equation 

with this heat flux. 

The wall boundary conditions for < > are based on the fact that at the wall, 

< u0 >=< v9 >= 0.0 

Therefore, the value of < vO > at the wall adjacent node is set by interpolation 
between the wall and the node next to the wall adjacent node. The value of < u9 > 
at the wall adjacent node is set equal to minus twice the value of < v9 > [39]. 

Since the velocity-temperature products < U{Uj9 > fall to zero at the wall, the 
near-wall values, hence, were made very small such that zero wall values is obtained. 
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Outflow Conditions 

Zero streamwise gradient conditions are used as the outflow conditions for all 
variables with the exception of the U momentum and temperature. 

The law of conservation of mass is used as the outflow condition for U momentum. 
The outflow temperature distribution is prescribed in accordance with Kays and 
Crawford [47] as 

dT = 2g" 

dx pCpU D c 


2.5.2 Turbulent Jet Flow 


Inlet Conditions 

In the self-preserved region of the jets, the U-velocity is assumed to follow the 
cosine curve. Therefore, a cosine function is used to describe the inlet U-velocity with 
the amplitude of the free stream velocity. V-velocity is assumed to be zero at the 
nozzle. 

The turbulent kinetic energy and turbulent energy dissipation is prescribed by the 
following relations: 



where l m is the minimum of the wall mixing length and the Nikuradse s formula. Both 
formulae are listed below with the wall mixing length relation first, and followed by 
Nikuradse’ formula. 


L = 0.41 y 

^ = 0.14- 
tn 


1 

to 

1 

O 

o 

05 


V r N J \ 

Tat/ 


( 66 ) 

(67) 


The Reynolds stresses are prescribed by using the Boussinesq viscosity correlations 


as 


— p < UiUj >= p t 


m duA 

dxj dxi ) 



The inlet temperature is given the value described in the experiment. 


( 68 ) 



Chapter 3 

Numerical Procedure 


3.1 Elliptic Code 

3.1.1 Discretization of Governing Equations 

The steady state, two dimensional governing equations of all dependent variables 
solved in the backward-facing step flows can be written in the following form 




d 

[r*l 

d 

-f" — — 

r d<f> ' 

dx 

dx 

dy 

dy 


+ S<j> 


(69) 


where <f> is the dependent variable £/, V, fc, e, < U{Uj > , .... This equation is divided 
into three general parts, namely convection of <f> , diffusion of (f> and source term of <f> 


Equation ( 69 ) is discretized to a linear algebraic equation before solution 
of dependent variables can be obtained. Before proceeding further, a few definitions 
of grid system must be made. Figure 2 outlines detail definitions of geometrical 
dimensions and locations. Letters P,E,W,N,S are the node P inside it’s control volume 
and it’s neighboring nodes East, West, North, and South, respectively. The letters 
e,w,n,s are the control surfaces of P control volume at east, west, north and south, 
respectively. 

Discretization of equation ( 69 ) is carried out by using the control volume ap- 
proach of Patankar [48, 49]. The final form of the discretization equation using central 
differencing scheme is as follow: 
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Figure 2: Definitions of geometrical dimensions and locations. 


ap4>p — ap<j)p aw<t>w + a,N<f>N + + b 


where 


( 70 ) 


= 

1*, A 1 t r a 

A^-Ax--,KA* 

a s = 

a>"^a* 

a E — 

^V-\pV,Ay 

a\v = 

^2™ a 1 

S^Aj, - -^Ay 

Cl p = 

a E + aw + a N + a s 


b = 


SuAxAy 


Although central differencing scheme has second order accuracy, it is not used 
in the computations because some of the coefficients, a,- might be negative at cer- 
tain iterations. This is violating the four basic rules which can result in physically 
unrealistic solutions. As a result, a new scheme called upwind differencing scheme 



scheme 

Formula for A(| P |) 

Central Difference 

Upwind 

Hybrid 

Power law 

Exponential 

1 - 0.5 | P | 

1 

|| 0, 1 — 0.5 | P | || 

|| 0, (1 — 0.5 | P |) 5 || 
1 P 1 / ( ex P 1 P 1 - 1) 


Table 3: A(| P |) for various differencing schemes. 

was introduced. This scheme has a set back because it only has first order accuracy. 
Spalding [50] proposed a scheme that combines the advantages of central differenc- 
ing and upwind differencing schemes and called it hybrid differencing scheme. With 
minor modifications, the general discretization equation can be written as: 


ap(j)p = a,E<j>E + aw<f>w + aN<f>N + as<f>s + b 


(71) 


where 


ajv = D n A(\P n \)+ || — P„, 0 || 
a s = D.A(\P,\)+\\F.,0\\ 
a E = D e A(\P e \)+ || — P e , 0 || 
a\v = D W A(\P W \)+ || F w , 0 || 
ap — g-e + aw T ojv + as — SpAxAy 
b = SuAxAy 


and 


r a 


D e = -r—Ay, D w = -^Ay, D n = ^ Ax, D s = 


2 n 


Ax^ 

P. = ^Ay, 

D e J ’ 


Ax 


A y n 


r 2s 

Ay 3 


Ax 


p.A-Ak Ay , 


P, = QA Ax 

F e = (p[/) e Ay, F w = (pU) w Ay, F n = (pV) n Ax, F s = ( P V) S Ax 


D, t 


D n 


( 72 ) 


the formulae for A(\P\) are listed in table 3 for various schemes. 
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Figure 3: Function A(| P | ) for three different schemes. 

Smce the goodness of different schemes are based on profile fitting with the one 
dimensional analytic solution, Amano [51] proposed a fourth order differencing scheme 
by expanding the exponential scheme of Spalding and found that this scheme gives 
better agreement with the exponential scheme than the hybrid scheme. Figure 3 
shows the graph of A(\P\) and the fourth order scheme is indeed better than the 
hybrid scheme. The formula of A(|P|) for the fourth order scheme can be written as: 


A( I P D =I1 °* 1 " b Pl + I2 |P|2 ~ ^l 4 II (73) 

3.1.2 Grid System 

The grid system used in this study is the so-called staggered grid system in 
which all scalar quantities are associated with every grid node (i.e. points where grid 
lines intersect), while the U-cell is displaced half a node leftward, V-cell half a node 
downward and < „„ >-cell both half a node leftward and downward. Figure 4 shows 
the grid, storage locations and control volumes of all dependent variables. The grid 
system is advantageous in solving the velocity field because the pressure gradients 



Figure 4: Staggered grid system and the grid storage locations used in elliptic com- 
putations. 

can be evaluated easily and velocities are conveniently located for the evaluation of 
convective fluxes [20]. 

3.1.3 Solution of Discretization Equation 

After the formulation of the discretization equation, a line by line iterative method 
is used to obtain solutions for all' dependent variables. The solution domain is given 
some initially guessed values which are improved upon from one line to the other. In 
this procedure, the values of 4 > on neighboring lines are assumed to be temporarily 
known. This, then reduces the unknowns to be solve to three and turns the discretiza- 
tion equation into a convenient algorithm sometimes called the Thomas algorithm 
or the TDM A (Tri-Diagonal-Matrix Algorithm). This algorithm is used along the 
North-South line and swept from west to east of the solution domain. A more detail 
explanation of the algorithm is given by Amano [52]. 



3.1.4 Pressure and Velocity Corrections 

At this stage, all dependent variables can be evaluated by the TDMA algorithm 
except the pressure field. The pressure gradient forms part of the source term in the 
momentum equations. If the correct pressure distributions of the solution domain are 
known, there will be no difficulty in the solution of the momentum equations. 

As mentioned before, the solution domain was given some guess values which most 
likely were not the correct distributions for most variables. This made the solution of 
the momentum equations impossible (due to incorrect pressure field). To correct this 
problem, the pressure field is indirectly specified via the continuity equation. When 
the correct pressure field is substituted into the momentum equations, the resulting 
velocities fields satisfies the continuity equation also. 

This indirect substitution of the pressure field is called the SIMPLE algorithm 
(Semi-Implicit Method for Pressure-Linked Equation). Patankar [48, 49] outlined the 
SIMPLE algorithm in its detail. 


3.2 Parabolic Code 

3.2.1 Discretization of Governing Equations 

The steady state, two-dimensional parabolic flows governing equations can be 
written in the general form of: 




dy\ 


+ s 0 


(74) 


where <f> is again the dependent variable (/, V, < U{Uj > ... The only difference 
between equation ( 74 ) and equation ( 69 ) is the omission of axial diffusion in 
parabolic flows. 

Equation ( 74 ) is discretized to a linear algebraic equation following Patankar’s 
finite-volume approach with hybrid differencing scheme for advection and central 
differencing scheme for diffusion terms, respectively. 

The final form of the discretization equation using hybrid differencing is as follow 


a P <f)p = ap(f>E + a\v4>w + aN<f>N + &s<f>s + b 


(75) 
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where 


a N = D n A(\P n \)+ || — F„,0 || 
g 5 = D„A(\P a \)+ || F s ,0 || 

a E — 0*0 

aw = || F wy Q || 

ap — ctp + aw + u;v 4“ &s — Sj^AxAy 

b — SVAxAy 

The definitions of all variables are given in equation ( 72 ) and A(|P|) is the hybrid 
differencing scheme given in table 3. The coefficient a# is set to zero because there 
are no influence from downstream nodes for parabolic flows. 

3.2*2 Grid System 

In this code, the grid system used is slightly different from the one used in elliptical 
computations. In this system, all scalar quantities including the Reynolds shear stress, 
< ur >, are associated with the grid node, while U-cell and V-cell are displaced half a 
node leftward and downward respectively. Figure 5 shows the grid storage locations 
and control volumes of all dependent variables. 

3.2.3 Solution of Discretization Equation 

After formulating the discretization equation, a line by line noniterative procedure 
is used to obtain solution for all dependent variables based on the known upstream 
values. A noniterative procedure is used as compared to the iterative procedure in 
elliptic computations because there are no downstream influence in parabolic flows. 


3.3 Boundary Conditions 

In this study, the grids are arranged such that cell boundaries coincide with 
edges of the solution domain. Figure 6 depicts a typical grid and solution domain 
ari'angement. Since some nodes lie outside of the solution domain, the solution is 
not influenced by these nodes. As a result, the coefficients for these nodes must be 
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Figure 5: Grid storage locations and control volumes used in parabolic computations. 


set to zero. In figure 7, point W lies outside of the solution domain. Therefore, 
at point P, the coefficient a.\y must be set to zero. At this point, the correct value 
may be added through the source term. Generally, there are three types of boundary 
conditions, these are the Dirichlet condition (Prescribed values of <j> B ), Neumann 
condition (prescribed flux at boundary) and Robin condition (boundary flux specified 
via a coefficient and condition of surrounding fluid). Apart from these, there is a so- 
called internal condition, where any internal node can be set to a desired value of <j>p. 
Amano [52, 53] outlined the mathematical formulations for the treatments of these 
conditions. 


3.4 Grid Generation 

An algebraic grid generation technique was developed to generate boundary 
conforming grids for two-dimensional duct flows. In this technique the grid point 
locations on the top and bottom boundaries are given, and the grid points in the 
interior of the domain are computed along straight lines connecting the corresponding 
top and bottom boundary grid points. The method described here allows for the 


Figure 6: Typical 





M 
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Figure 7: 
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generation of grids with the grid line spacing next to the walls kept at a constant user 
specified value. 

If (TV _ 2) grid points are to be distributed along a straight line between the grid 
points (xi, yi) and (x N , y N ), and if the spacing between the points (*i, yi) and ( x 2 , y 2 ) 
and the spacing between the points (x^-i, Vn-i) an( l ( x Ni Vn) is to a user specified 
value, (As)i, then the expansion factors in the x and y directions are given as: 


Px = 



AX(fc ~ 1) 
2(Ax)i 


2 

77=T 


(76) 


Py = 



A Y{Py - 1) 
2(Ay)i . 


2 


where 


AX = xjv — xi 


(77) 


(78) 


Ay = vn - 2/1 


(79) 


(Ay)i = sign(AY) 



(80) 


(Ax)i = (Ay) 


AX 

1 AY 


(81) 


Since the equations for the expansion factors p x and p y cannot be solved explicitly, it 
is necessary to iterate. Fixed point iteration using the expressions in the form shown 
above was found to converge within a few iterations, and was used to generate all the 
grids shown here. Once the expansion factors are known, the grid point locations are 
given by: 


x n 


== X\ + (Ax)i 


( lx)"' 1 ~ 1 
Px~ 1 


for n < — ^ — (82) 


yn ^i d - (Ay)i 


(Ml ~ 1 

A-i 


for n < 


N+l 


2 


(83) 
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(8 \ N ~ n — 1 N A- 1 

x n = x N - (Ax) r ^ — for n > — - — (84) 

(8 ) N ~ n — 1 N 4 - 1 

Vn = VN - (±yh —p ~1 forn> — —— (85) 

The method shown here was found to be very efficient, and could be used for many 
different geometries. 

3.5 Coordinate Transformation 

Equations ( 69 ) are transformed from the Cartesian coordinates (x, y ) into general- 
ized curvilinear coordinates (£, rj). Consider the general steady state two-dimensional 
transport equation in the following form: 

{pU(j>)x + (pv<t>) y = (r<k)* + (!>„)„ + R^{x,y) (86) 

where T represents a diffusion coefficient and R $ shows a source term of the transport 
equation of a dependent variable <j>. Equation ( 86 ) can be written as 

j[W)f+W),] = J - W,)] +7[j(^-W<)] +S.K,i»I87) 

where S^^rj) is the transformed version of the source term i?^(x,y) and where 


J = - x v yt (88) 

a = x l + yv ( 89 ) 

P = xt x v + y$y v ( 90 ) 

7 = x| + y£ (91) 

In Eq.( 87) the contravariant velocities U and V are given as: 

U = Uy v — Vx n (92) 

V = Vxf. - Uy £ ( 93 ) 


Equation ( 87 ) is discretized by using the control volume method described in the 
following section. 



3.5.1 Numerical Procedure for Transformed Coordinate Sys- 
tem 

Formulation and discretization of all transport equations were performed by using 
the conventional control- volume approach of Patankar [48], by breaking each equation 
into diffusion, convection, and source terms. The system of equations were made 
linear so that they could be solved iteratively by the tridiagonal matrix algorithm. 
After several numerical tests, it was observed that the variation in results of mean 
velocity profiles lies within 2% when the grid is changed from 52 x 52 to 62 x 62 for 
the computations of an angles backward facing step flow (Figure 8). Therefore the 
grid independent state was assumed to be attained with the grid of 62 x 62. 

3.5.2 Pressure Correction Algorithm 

Three pressure correction algorithm were tested: SIMPLE [55], SIMPLEC [56], and 
PISO [57]. Here we review the approximations made by the three algorithms. 

SIMPLE 

The SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm 
can be derived by first obtaining a discretized form of a velocity correction equation. 
The U-component velocity correction equations: 

A e U' v = EA„ t UU - + ytP', (94) 

where A’s are the influence coefficients, and U’ and P’ are the velocity and pressure 
corrections. A with the subscript nb denotes the neighbor coefficients. If the above 
equation were used in the derivation of the pressure correction equation, an unman- 
ageable equation would result. Instead, it is assumed that the velocity correction 
at a point is not affected by the velocity corrections of its neighbors, and thus the 
summation term in equation ( 94 ) is neglected. 

SIMPLEC 

The SIMPLEC (SIMPLE Consistent) algorithm attempt to use a more consistent 
approximation, based on the magnitude of the terms in the velocity correction equa- 
tion. Instead of neglecting the velocity corrections at the neighboring points, the term 
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T,A nb Up is subtracted from both sides of equation ( 94 ). This leads to the following: 

(A p - ZA nb )U' p = ZA nb (U^ b - %) - y v P[ + y (95) 

Again, if the above equation were used to form the pressure correction equation, an 
unmanageable equation would result. In the SIMPLEC algorithm the summation 
term on the right hand side is neglected. 

PISQ 

In the PISO (Pressure Implicit with Splitting of Operators) algorithm the pressure 
and velocities are corrected by a series of steps. In the incompressible form used here, 
it is assumed that two pressure and velocity corrections are sufficient. The first 
pressure correction equation of PISO is identical to that used by SIMPLE. A second 
corrector step is done to ensure that the continuity and momentum equations are 
satisfied at the end of each iteration. The second corrector step requires the solution 
of a second pressure correction equation, and thus the computation time per iteration 
will be longer for PISO than for either of the other algorithms presented here. 

Due to the grid being non-orthogonal, the pressure correction equations contain 
cross derivatives, which lead to a nine-point formulation (ie. the pressure correction 
at a point P is a function of the pressure corrections at its neighbors, N, S, E, W, 
NE, NW, SE, and SW). For this reason it is necessary to either neglect the cross 
derivatives, incorporate them into the source term, or use a nine-point solver, since 
Shyy et al. [57] found that incorporating the cross derivative terms in the source 
term had no advantages over the method of neglecting the cross derivative terms. 
It is also reported by Ando et al. [58] that the use of a nine-point solver ensured 
stability and robustness of the computational method. For the reasons mentioned 
above, the semi-implicit solver of Peric [59] was employed in this study to solve the 
pressure correction equation. 


3.5.3 Velocity Correction Algorithm 

The Cartesian velocities U and V were corrected and contravariant velocities U and V 
were calculated using the new values of U and V rather than computing the Cartesian 



velocities and correcting contravariant velocities. This is because it was found that 
the former method gives more accurate results and is more efficient than the latter. 
The velocity correction equations for U and V are: 


u = u* + U‘ 

(96) 

v = V* + V' 

(97) 

Substituting the velocity corrections U' and V': 

U = U* + B% + C% 

(98) 

v = v* + b v p [ + 

(99) 

where 

B u * B v “ Xt} * C u — ^ * C v — X ^ 

~ Al' ” A v ' A u 

(100) 


and where A't and A v v are the coefficients of U p and V p , respectively. U and V were 
calculated using the definition of the contravariant velocities. 



Chapter 4 


Results and 


Discussion 


This chapter presents the numerical results obtained from the modeling of turbu- 
lence phenomena in backward-facing step and jet flows. 

The experimental data of Driver and Seegmiller [2] and Chandrsuda and Bradshaw 
[5] were used to validate the turbulence models developed in chapter 2 for backward- 
facing step flows while, the experimental data of Antonia et al. [8, 9, 10, 11, 12, 13] 
, Heskestad [14] and Gutmark and Wygnanski [16] were used to validate the model 
for jet flows computations. 

The iterative solution procedure is terminated when the maximum value of the 
relative residual sources of U, V and mass balance falls below 1%. However, the 
computations of the triple products are terminated when the relative residual sources 
fall below 3 x 10“ 8 for < u x UjUk > and 5 x 10" 9 for < U{Uj9 >. 

The complete process of solving the momentum, temperature and related turbu- 
lence products equations takes approximately 60 minutes of CPU time on a UNIVAC 
1100 computer. The CPU time varies with Reynolds number of the flow. 

Figures 9 and 10 show numerical grids used in the computations of backward- 
facing step flow. The figures shown here are the actual computation domains used for 
experiments of Driver and Seegmiller and Chandrsuda and Bradshaw, respectively. 
These are 62 x 62 mesh grid system with the grid expanding linearly at the rate of 
2% and 3% in the axial and transverse directions, respectively. The height of the 
channels are magnified five times for better visualization. 

Figures 11 and 12 show the velocity variations inside the solution domain for the 
experiment of Driver and Seegmiller and Chandrsuda and Bradshaw, respectively. 
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Figure 9: Grid system used in the computations of backward-facing step flow for the 
experiment of Driver and Seegmiller (5:1 magnification radially). 



Figure 10: Grid system used in the computations of backward-facing step flow for the 
experiment of Chandrsuda and Bradshaw (5:1 magnification radially). 
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Figure i-F Velocity variations inside the solution domain for the experiment of Driver 
and Seegmiller. 



Figure 12: Velocity variations inside the solution domain for the experiment of Chan 
drsuda and Bradshaw. 
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Figures 13 and 14 show variations of the computed Reynolds stresses for the 
solution domain. 

Figures 15 and 16 show the Low-Reynolds and high- Reynolds number model’s 
third-moment variations, respectively. It can be seen that the high-Reynolds number 
model predicted higher levels at the near-wall region. 

The solutions are obtained by using a 75 constant grid mesh system across 
the traversing direction of the jet and computations are terminated when the desired 
axial location is solved. 

The complete process of solving the momentum, Reynolds stresses and tempera- 
ture equations takes about 5 minutes of CPU time on a UNIVAC 1100 computer. 

Figure 17 shows the comparisons of computational results with the measured 
velocity profiles of Heskestad [14], Gutmark and Wygnanski [16] and Antonia et 
al. [9] at self-preserving region. Figure 18 shows the comparison of the computed 
temperature profile with the data of Antonia et al. [9] - The agreements are very 
satisfactory in all three cases. 

Results by Using Different Pressure-Strain Correlation Model 

Figures 19 and 20 show the comparisons of computed Reynolds stresses with 
the experimental data of Heskestad [14] and Gutmark and Wygnanski [16]. Three 
different models are used, and the transport equations model with the pressure strain 
correlations developed here (see sec 2.3). The new model gives excellent predictions 
for the Reynolds stresses < uu > and < vv >. However, this model always give 
higher turbulent shear stress, < uv >. 

Comparisons of elliptic flows are shown in Figures 21 to 27. Figure 21 shows 
the computed mean velocity profiles and are compared with the experimental data of 
Chandrsuda and Bradshaw and of Driver and Seegmiller, respectively. 

Figures 22 - 27 show results of the Reynolds stresses for the same cases as men- 
tioned above. It is obseved that the new model considerably improves the predictions 
for all the components of the Reynolds stresses except < vv >. The component of 
< vv > are underpredicted by employing the present model when compared with 
other models. However, it is also noticed that the experimental data for < vv > are 
much higher in the cases of backward-facing step than those for jet flows. Thus, the 
predictions by the new model may be more consistent with the data of jet flows than 
those of step flows. 




Figure 13: Reynolds stresses 
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Figure 18: Temperature profile - Comparison of the computed temperature profile 
with the data of Antonia et ah. 

Flows Through Irregular Boundary Ducts 

Calculations were performed for laminar flow (Re = 500) in a kidney shaped chan- 
nel with a 25x32 grid using the SIMPLE, SIMPLEC, and PISO pressure correction 
algorithms. The calculations were considered to be converged when the normalized 
mass and momentum residuals were reduced to less than 0.01 %. The computed veloc- 
ity vectors for the flows through the kidney shaped channel and through a diaphragm 
pump chamber are shown in Figure 21. 

Computed velocity vectors through a 180° turn around duct is shown in Figure 

22 . 

A comparison of the three pressure correction algorithms is shown in Figure 23. It 
is shown that both SIMPLEC and PISO require approximately 45 % fewer iterations 
than SIMPLE to meet the above mentioned convergence criteria. Since one PISO 
iteration takes longer to complete than either one SIMPLE or SIMPLEC iteration, 
comparisons were also made on the basis of work units (WU), where we define one 
WU to be equal to the time required to complete one SIMPLE or SIMPLEC iteration. 
Figure 24 shows the comparison of the three algorithms on the basis of work units. 
Here it is seen that PISO requires 29 % less computation time, and SIMPLEC requires 
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Figure 21: Mean velocity profiles downstream from step. 
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Figure 23: Reynodls-stress profiles downstream from step. 
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step. 


Figure 25: Reynodls-stress profiles downstream from 
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Figure 26: Reynodls-stress profiles downstream from step. 
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Figure 28 : Velocity vectors in a kidney shaped channel (top) and in a diaphragm 
pump chamber (bottom) 





Figure 30: Comparison of three different pressure correction algorithms. 


41 % less computation time than SIMPLE. 




normalized mass residual 



WORK UNITS 


Figure 31: Comparison of three different pressure correction algorithms. 




Chapter 5 
Conclusions 


After completing the studies presented in the previous chapters, a few observations 
is made and presented below. 

There are a few advantages in using the full transport equations to solve 
turbulent flows. Presented below is a list of these advantages. 

1. The Reynolds-Stress model developed in this study predicts the separating and 
reattaching shear flows properly. 

2. The transport equations for the third-moment of turbulence predicts the rapid 

changes of third-order turbulence velocity fluctuating tensor in the reattaching 
and recirculating flow regions. 

3. The transport equations are more superior to the algebraic equations in the 
predictions of turbulence quantities because the convective and diffusive effects 
neglected by the algebraic models are accounted for by the transport equations. 

4. The Low- Reynolds number model of third-moment of turbulence which pro- 
motes the dissipation effects of the third-moment in the near-wall region, im- 
proves the predictions of third-moment considerably and gives more universal 
results than the algebraic equations. 

5. The newly developed parabolic code predicts the turbulent jet flows well, and 

further investigations and testings should be carried out to further test the 
abilities of this code. 



6. The model predicted the mean velocity profiles to within 5% of the measured 
values. 

7. The predictions of Reynolds stresses is very sensitve to the choice of pressure- 
strain correlation. Presently, the model of Amano et al. [35] produces the best 
predictions. 

8. This code is more preferable than the code of Patankar and Spalding [31] because 
no coordinate transformations is needed. 

9. The SIMPLEC and PISO algorithms require far fewer iterations than the SIM- 
PLE algorithm. This shows that the Consistent” approximation used in SIM- 
PLEC and the operator splitting approximation used in PISO are better than 
the assumption used in SIMPLE, namely that the velocity correction at a point 
is not affected by the velocity corrections of its neighbors. 

10. The SIMPLEC algorithm proved to be more efficient than either PISO or SIM- 
PLE for the case studied. 

11. The PISO algorithm uses more memory than either SIMPLE or SIMPLEC, 
but has an advantage over the other algorithms in that it can be used for non- 
iterative time-dependent solutions. 

12. The increased efficiency of SIMPLEC and PISO compared to SIMPLE would 
probably be more dramatic if a finer grid had been used, since the performance 
of SIMPLE degrades dramatically as the grid is refined [61]. 

Finally, in order to predict the turbulence quantities accurately, the solution of 
the full transport equations of each individual quantity is needed. 
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